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Abstract 

Assuming that nuclear matter can be treated as a perfect fluid, we study the propagation of 
perturbations in the baryon density. The equation of state is derived from a relativistic mean 
field model, which is a variant of the non-linear Walecka model. The expansion of the Euler 
and continuity equations of relativistic hydrodynamics around equilibrium configurations leads 
to differential equations for the density perturbation. We solve them numerically for linear and 
spherical perturbations and follow the propagation of the initial pulses. For linear perturbations 
we find single soliton solutions and solutions with one or more solitons followed by "radiation". 
Depending on the equation of state a strong damping may occur. Spherical perturbations are 
strongly damped and almost do not propagate. We study these equations for matter at finite 
temperature. 
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I. INTRODUCTION 



Over the last decades hydrodynamics of strongly interacting systems l|, |2|, y, |J] has been 
applied to cold nuclear physics, to low and high energy nuclear reactions and to phenonema 
taking place in dense stars. Recently hydrodynamical models became more sophysticated 
and received more support from experimental data, in particular from the measurement of 
elliptic flow at RHIC [5]. Whereas other approaches can give satisfactory descriptions of the 
measured transverse momentum distributions, when it comes to elliptic flow, there are not 
many options other than hydrodynamical models. There is compelling evidence that we have 
seen "the perefect fluid" at RHIC. This evidence might be signiflcantly reinforced by the 
observation of waves. Waves in a hadronic medium are produced in many physical situations. 
In fact, there are already some indications that these waves have been formed. In relativistic 
heavy ion collisions we may have hard parton - parton collisions in which the outcoming 
partons have to traverse the surrounding fluid to escape and form jets. Their passage may 
form Mach shock waves {g!], which will affect the transverse momentum distribution of the 
observed flnal particles. These "Mach cones" may have been observed at RHIC 7|, |8j. 

Under certain conditions waves may form solitons. Therefore we can go a step further 
and look for solitons in a hadronic medium. In the RHIC scenario, for example, the same 
supersonic motion that generates conical shock waves may also generate solitons. Whether 
or not this happens, depends on details of the equation of state and on the approximations 
used in the hydrodynamical description of the motion. Another scenario where solitons may 
appear is in the core of dense stars. Here perturbations in the baryon density may be caused, 
for example, by interactions of neutrinos with the baryonic matter. In a pioneering series of 
works on soliton formation in nuclear matter o], [lo|, 11 1 it was suggested that in nucleon - 
nucleus collision at intermediate bombarding energies (~ 50 — 200 MeV) the nucleon may 
be absorbed by the nucleus (treated as a fluid at rest) and propagate as a localized pulse of 
baryon density. 

In this work we study the propagation of sound waves in dense and hot hadronic mat- 
ter. More speciflcally we consider the propagation of perturbations in the baryon density. 
These perturbations may generate ordinary waves, shock waves and also Korteweg - de Vries 
(KdV) solitons. Starting from the equations of relativistic hydrodynamics at zero and flnite 
temperature and in Cartesian (x, t) and spherical (r, 6, (p, t) coordinates, we derive differen- 
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tial equations and find their numerical solutions. The equation of state is derived from a 



relativistic mean field model of the Walecka type [12|, [131, . We discuss the features of the 
solutions and the role played by the microscopic interactions in the shape and propagation 
of the sound waves. 

In a previous work 15| we have studied the formation and propagation of KdV solitons 
in cold nuclear matter. We found that these solitary waves can indeed exist in the nuclear 
medium, provided that derivative couplings between the nucleon and the vector field are 
included in the interaction Lagrangian. For this class of equation of state (EOS), which 
is quite general, perturbations on the nuclear density can propagate as a pulse without 
dissipation. 

During the analysis of several realistic nuclear equations of state, we realized that, very 
often the speed of sound Cg is in the range 0.15 — 0.25. Compared to the speed of light these 
values are not large but not very small either. This suggests that, even for slowly mo ving 
nuclear matter, relativistic effects might be sizeable. We investigated these effects in 16 | 
and in 17l |. 

The propagation of density pulses might be relevant for the astrophysics of dense stars. 
This motivated us to extend our results to the spherical geometry. In p^] combining the 
Euler and continuity equations in relativistic hydrodinamics in spherical coordinates, we 
have obtained for the first time an equation similar to the KdV equation. Spherical KdV 
- like equations have been found before in other contexts, as, for example, in the study of 
nonlinear waves in dusty plasmas [19|, |20| . 

In the present work we reexamine all our previous works, looking for the numerical 
solutions of the previously encountered differential equations. In the linear case we compare 
the analytical solution with the numerical one and study the sensitivity of the solution to 
the initial conditions. In the spherical case there is no analytical solution. Our numerical 
solution could be compared to the one found in [lol, Q . We have extended the formalism to 
finite temperature and, with the new equation of state, derived differential equations which 
are temperature dependent. We have also studied the limiting case where the differential 
equations generate shock waves. 

The text is organized as folows. In section II, for convenience, we collect some useful 
formulas for hydrodyamics. In section III we present the equation of state obtained with 
our model. In sections IV we derive the spherical KdV-like equations and in section V we 
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present and discuss the numerical solutions of these equations. Finally, in section VI we 
present some conclusions. 



II. RELATIVISTIC HYDROYNAMICS 

In this section we review the main equations of relativistic hydrodynamics. In natural 
units (c = 1) the velocity four vector m'" is defined as: 



u 



(1) 



where 7 is the Lorentz contraction factor given by 7 = (1 — f^)^^/^. The velocity field of 
matter is v = v{t,x,y,z) and thus u'^u^, = 1. The energy-momentum tensor is, as usual, 
given by: 

T^,y = + p)ut,Uy - pg^,y (2) 

where e and p are the energy density and pressure respectively {qqq = —gu = 1 and g^y = 
if /i 7^ z/). Energy- momentum conservation is ensured by: 



= 



(3) 



The projection of ([3 
of Euler equation [21 



onto a direction perpendicular to u'^ gives us the relativistic version 
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23|: 



The continuity equation for the baryon number is 2l|, |22, 1231]: 



(4) 



d.jB" = 



(5) 



Since Jb'^ = Ps, where pb is the baryon density, the above equation reads 

d 

-QlipBl) + V • {pBlv) = 



(6) 



or 



dpB 2 I 9v _^ -L 
^ + 7^P^^ + ^-V. 



V • {PBV) = 



The enthalpy per nucleon is given by 22|: 



(7) 



dh = Tds + Vdp 
4 



(8) 



where V = 1/ pb is the specific volume. T and s are the temperature and entropy 
density respectively. For a perfect fluid {ds = 0) the equation above becomes dp = psdh 
and consequently: 



The Gibbs relation is 24l |: 



^ dp dh 



e + p = PbPb + Ts 

where /ib is the baryochemical potential. Inserting ([9]) and ffTOj) in (j4]) we find: 

PB 



\- (v ■ V)v = 



The enthalpy per nucleon can also be calculated with the expression 

dE 



Ul: 



h = E + pb 



dp 



B 



where E is the energy per nucleon given by: 



(9) 



(10) 



(11) 



(12) 



E = — 

PB 



which, inserted into f[T2|) yields 



h 



de 



dp 



B 



(13) 



It is clear that the "force" on the right hand side of (fTTI) will be ultimately determined by 
the equation of state, i.e., by the function e^pb)- 



III. EQUATION OF STATE 



Equation (jTTj) contains the gradient of the derivative of the energy density. If e contains 
a Laplacian of p^, i.e., e oc ... + ...V'^pb + then ( fTTj) will have a cubic derivative with 
respect to the space coordinate, which will give rise to the Korteweg-de Vries equation for 
the baryon density. The most popular relativistic mean field models do not have higher 
derivative terms and, even if they have at the start, these terms are usually neglected during 
the calculations. 

As in [isl we shall use a variant of the non-linear Walecka model [3] given by: 

c = Cqhd + ^i>{d,d''v^)r^ (14) 
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with 



bcf)-^ C(f)^ 1 



4 A'""" ■ 2- 

where F^u = d^Vu—duVfj,. As usual, the degrees of freedom are the baryon field ip, the neutral 
scalar meson field and the neutral vector meson field V^, with the respective couplings 
and masses. The second and new term in f[T^ is designed to be small in comparison with 
the main baryon - vector meson interaction term g^ip'^^jy^il). Because of the derivatives, it 
is of the order of: 

ZL^J^^0.12 (15) 

my my 

where the Fermi momentum is kp — 0.28 GeV and my ~ 0.8 GeV. The form chosen for 
the new interaction term is not dictated by any symmetry argument, has no other deep 
justification and is just one possible interaction term among many others. It is used here 
as a prototype. At this stage our main interest is to explore the effects of these higher 
derivative terms, which may generate more complex wave equations. The parameter d is 
free and will act as a "marker". Setting d equal to zero will switch off the new term. On 
the other hand d = 1 means that the coupling gy is the standard one. Other values imply 
a correction in this coupling. 

The mean field approximation means that ^< >= S^qVq and (f) — (f) >= 
00- With the above Lagrangian and the corresponding Hamiltonian we can, following the 
standard procedure, write the partition function of the system and calculate the energy 
density, the pressure and entropy density, which, for symmetric nuclear matter, are given 



by |25|: 



p = 7r—2pB' - ?^(^ - Mr - - Mr - T^(M - Mr 

2my^ 2gs^ Sgs-^ 4:gs* 

-Tj^s I !^ln[{l-n^,{T,u))]+ln[{l-n^,iT,u))] 

-PbV pb + —^PB (17) 



and 



7s 



d^k{nAT,i')ln nt{T,ij) + l-nAT^v) In l-nAT^v) 



+ 



-\-n^{T,v)ln nj:{T,v 



+ 



n^(T, v) In I- n^{T, v) 



(18) 



where 7^ = 4 is the degeneracy factor and M* is the effective nucleon mass ( M* = M—gs4>o) 
given by: 



M* = M - 



9s Is 
ms^ (27r)3 



M 



+ 



9s^ 



and 



with 



PB 



Is 



9s^ 



(27r)3 
n^{T, u) = 



1 



1 + (,(h+~u)/T 
1 



(19) 
(20) 

(21) 

(22) 
(23) 
(24) 



1 + (,ih++^)/T 

V = - 9vVo 

The last two terms of (fT6!) come from the new interaction term in (fT4|) . Baryon number 
propagation in nuclear matter has been studied in 

dpB 



dt 



261] with the help of the diffusion equation: 
D V^pB (25) 



where the diffusion constant D was numerically evaluated as a function of density and 
temperature and found to be D ~ 0.35/m at densities comparable to the equilibrium nuclear 
density and temperatures of the order of 80 MeV. This number is small compared to any 
nuclear size scale and can be interpreted as indicating that density gradients do not disappear 
very rapidly in nucler matter. Using fl25l) twice in the last term of (fT6l) it can be rewritten 
as: 



dgv ( d'^pB 
XPb' 



my 



dgv d I 2 
-ps— DV pb 



'dt 



dgv" 



PbD'[V\V'pb)] (26) 



which, in the context of the present calculation, can be neglected because V^(V^pb) << 
(V^ps). With this last approximation, the final form of the energy density is given by ( fT6|l 



without the last term. Of course, the same argument holds for the pressure, which will be 
given by f|T71) without the last term. 

When T = the energy density reduces to: 

.2. 2gs^ 



2m 



V 



+ — — (M-M*)^+ , 



(27) 



From f|T6l) + f|T7|) we can check that the Gibbs relation ffTOj) is fulfilled. The speed of sound 
Cs is given by: 



9p 



(28) 



The numerical inputs for the above formulas for T = were taken from and are 
shown in Table I. In the Table the incompressibility K, the effective nucleon mass M*, the 
speed of sound Cg and the saturation density po are calculated. For T > the sign of the 
parameter c is reversed. In order to test our routines we have reproduced the results shown 
in Table I, but in what follows all the results will be obtained with the parameter set NLl. 





QHD 


NLl 


NL3 


A^L3 - // 


NL - SH 


K(MeV) 


545 


211 


272 


272 


355 


M(MeV) 


939 


938 


939 


939 


939 


msiMeV) 


500 


492 


508,2 


507,7 


526 


mviMeV) 


780 


783 


782,5 


781,9 


783 


9s 


8,7 


10,14 


10,22 


10,2 


10,4 


9v 


11,62 


13,28 


12,87 


12,8 


12,9 


M*/M 


0,56 


0,57 


0,6 


0,59 


0,6 


poifm-^) 


0,19 


0,15 


0,15 


0,15 


0,15 


b{fm~') 





+ 12,17 


+ 10,43 


+ 10,4 


+6,9 


c 





-36,26 


-28,9 


-28,9 


-15,8 


Cs 


0,25 


0,16 


0,18 


0,18 


0,2 



Table I: Numerical inputs [12] for the equation of state. K, M*, pQ and Cg are calculated. 
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IV. DERIVATION OF THE SPHERICAL KDV EQUATION 



A. Zero temperature 



This section contains the derivation of the spherical KdV equation. The general scheme 
is the same as the one used in 15| for the one dimensional Cartesian problem. However, as 



it wil be seen, there are new details, which deserve the discussion presented below. We shall 
be concerned only with problems which have spherical symmetry. Therefore the continuity 
equation ([7]) and Euler equation (fTTj) will have only radial components and the derivatives 
with respect to angles will vanish. 

Cold nuclear matter exhibits the saturation property, i.e., the energy per nucleon as a 
function of the baryon density {E = e/pb) has a minimum. Thus we start with (p7|) and 
impose the saturation condition: 

_d_ 

dps \PB 







(29) 



Pb=Pq 

We then perform a Taylor expansion of E around the equilibrium density po up to second 
order: 



E{pb) = E{po) + 



Ifd^E 



2\dp 



B 



{pB - Pof 



(30) 



As in 



Pb=Po 



the density pb and its gradient Vp^ are treated as independent variables. In- 



serting the above expression into (fT2l) and using the relation 



d^E 



Pb=Po 



we obtain: 



9v 

o 2P0 + 

zniv 



ms 
2po 



2 r 



(M* - M) 



9s 



+ 



Mcs 
Pi 



Is 

(27r)3po 



(31) 



{M -M*f + 



^9s^Po 



9v' 



^9s^Po 



(M - My + d^V'pB 



my" 



2po^ 



■(3pij^ -4pBpo + po^) 



(32) 



With the above expression we compute the derivatives: 



dh 
dr 



2Mc,2 dpB , 3Mc,2 dpB , , 9v^ d^PB 
H —PB— \- «- 



Po 



dr 



Po' 



dr rny^ dr^ 
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and also: 



dh 
'di 



jYiyi 7- Qf2 my^ dr 



2Mcs'^dpB 3Mcs'^ dpB .gv^dfd^pB 



(33) 



Po dt 



+ 



2 'PB- r,. 



+ d 



mv^ dt \ dr"^ 



+ 



my*^ r dt\ dr J 



Inserting fl33|) and fl34j) into f[TT|) we find: 



9f 9f 
St dr 



Pb 



fSMcs 



Po' 



-PB- 



2Mcs 
Po 



+ d 



2gv^ \d^PB IdpB 



rmv 



dr"^ 



r dr 



dpB ^ ^dpB 
dr dt 

, d ( dpB 

+ V 



+d 



dt V dr 



3X1 
f ^ dp 



d^pB , d ( d'^pB 



dr^ dt V dr"^ 



B 



r dr 



We now rewrite (|35|1 e ([7j) in terms of the dimensionless variables: 



(34) 



+ 



(35) 



. PB 
P = — 

Po 



V 



(36) 



where po is the equilibrium baryon density given in Table I and Cg is the speed of sound. 
The Euler equation becomes: 



dv 
dt 



Cs— + Cs V 



dv 



ic^H'^ - 1) r /3Mc 2 



dr pb 
d ( d'^p 



I- 



Po' 



-Pop 



2Mcs 
Po 



Po 



dp ^dp 
dr ^ dt 



+ d jPo 



niv 



d^p 
dr^ 



c.v- 



dt \ dr"^ 



2^y2 \d'^p _ 1 dp 

J, 



, -JV 

d jPo 

rmy* 



d f dp\ Cs V dp 
r dr 



and the continuity equation becomes: 



:i-c.,vi 



dp ^dv ^dp 2cspv 

— + Csp^ + CsV— H 

dt dr dr r 



dv 



dv 



+ CsPV[^ + CsV^^ 



(37) 



(38) 



We next define the "stretched coordinates" ^ and r as in |9l. llOl. Illl. l27l|: 



1/2 



[r - Cst) 



r = a 



3/2 



c.t 



(39) 



R ' R 

where i? is a size scale and a is a small {0 < a < 1) expansion parameter. The derivatives 
become the following operators: 

d d d^ a d^ d^ 



dr 



R di ' 9r2 R^de 



dr^ i?3 Q^Z 



and 



d _ a^/^Cs d a^d_ 
dt~ R di^ R dr 



(40) 



10 



From (139!) we can see that: 



and thus: 



R^a + Rt 



a 



3/2 



1 



r2 {R^a + Rrf 

In the i — T space ( 157|1 reads: 

f I f 1-^2 

R R dr ' 



(41) 
(42) 



1/2 dv 



(C^-O^ - 1) r /3Mc 



2Mc, 



, 1 2 "^0/^ ^0 



aV2c 2 



-V^ + zr^V-^ 1 



+d -po 



my" 



cr 



3/2 ^3^ 



V 1 V 

R R drjyR^de 



+ 



+ 



+d 



3/2 



mv'^iR^a + Rt) 



_ g^/^ dp 
R^dC (R^ct + Rt) R d^' 



+ CsV 



R R dTj \ R dij {R^a + Rt) R 

and the continuity equation (!38l) reads: 



(43) 



(1 - cg'v') - 



o^l'^Cs dp cr^/^Cs dp 



3/2, 



R di 



R dT 



+ Csp- 



1/2 



R di 



cr^/^ dp 2csVpo 



3/2 \ 



+ 



+ Cs pv 



o^l'^Cs dv cr^/^Cs dv 

H 7^ -TT- + CgV 



R dC R dT 
We then expand ( l36l) around the equihbrium values: 



R d^ (R^a + RT) 
a^/^dv\ ^ 



+ 



R 



p = 1 + api + a p2 -\ = 1 + pi + p2 + ... 



(44) 



(45) 
(46) 



With this expansion ( 1431) becomes: 



R 



f7-3/2 2 n 

+ (J t;2 + . . . ) H — (^^^i + ^ ^2 

R dT 



+Cs^——{(TVi + a'^V2 + ... )i:-{crvi + a'^V2 + ... 
R dt, 



[cs\av, + a^V2 + ...)-!] 



X 



P-B 

■po(l + api + a p2 + . . . ) - 



Po 



Po 



0" 



(l + Cipi + (7^2 + •••) + 
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+ 



R 



(cTf 1 + (T^f 2 



d 
d 



(1 + (7P1 +(7^2 + •••) + 



(crWi+cr^W2 + ---)7r(l + ^Pl+^V2 + 
it ar 



+ 



my* 



3/2 ^3 



(1 + api + (7^2 + ...) + 
a \ / a 



^(a., + a.2 + ...)^+ 



+d 



mv^iR^a + Rt) 



R^oe 



(1 + api + cr^p2 + • • • 



+ 



R^Oe 



(1 + api + a>2 + ...) + 



^3/2 ^1/2 ^ 



(1 + (Xpi + 0-^2 + ...) + 



(R^a + Rr) R 



+ Rt) 
and (jH]) becomes: 

[1 - Cs'^icrvi + a'^V2 + ■ ■ ■ ) 



R R dTj\ R di 

(1 + api + o-''p2 + . . . ) 



(1 + api + o-^p2 + . . . 



R 



(47) 



+ 



R dr 



R 

(1 + api + a^p2 + ...) + 



(l + (jpi + aV2 + •••) + 



0-1/2 ^ 

+Cs{l + api + 0-V2 + • • • + ^^^2 + •••) + 

+C5(o-Ui + a'^V2 + ■ ■ ■ )— + (^Pi + 0-^2 + ...) + 



2cs{avi + a^V2 + ...){l + api + a^p2 + ... )cr^^^\ 



+ 



{R^a + Rt) 

+Cs^(l + api + cr^p2 + . . .){avi + (T^f2 + . . . ) 
a3/2c, d 



+ 



a'/^Cs d 
R 



{avi + a'^V2 + ...) + 



{avi + a'^V2 + ...) + Cs{avi + a'^V2 + . . ■)—^^{(^Vi + a'^V2 + . . . 



R dr' ' ' / . .V . . z . ' R 

Since cr is small we go only up to second order. Therefore (H7j) and (HHl) turn into: 



(48) 



2Mcs 
PO 



-)Po 9pi 



Pb Cs 



+ a' 



dv2 dvi dvi 



'3Afc 



2Mcs 



— 2— Po 



■)P0 9p2 



PB Cs 



3Mcs'^ „ 2 



PB 



Po' 9pi (^ 
Pl- 



Po" 



Pb 



5^ ^ VPB C.2/?2 y 5^3 



Si' 
(49) 
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and 

■ dpi dvi\ 2 



dv2,dpi_ 2 dvi dp2 , 



9pi Stii 2 
+^1^ + /'i^ + 77 — ] — 7^1 







In the last term of the above expression, since < o" < 1, we shall assume that r > and 
make the approximation 

2 2 

= - (50) 

and hence: 

/ dpx dv{\ 2fdv2 ^ dpi 2 dvi dp2 dp, dvi 2 \ 

Since the coefficients in the above series are independent of each other we get a set of 
equations. From the terms proportional to a in (149|) and (15T|) we find: 

Pi = vi (52) 

and also 



and therefore 

/is = M (54) 

In fact, in fl52p we might have an integration constant. However, as it was shown in for 
the one dimensional Cartesian case, this would not change the results significantly. For our 
purposes it is enough to consider ( l52i) . keeping in mind that it is only a particular solution 
of the problem. From the terms proportional to cr^ in ( H9l) and (15T1) . with the help of (!52ll 
and (153|) . we find: 



2 

dpi ^ dpi ^^ dpi ^ ^2^ dpi fd^po\d^pi 

r3 



9pi 2 ^Pi , dpi dpi 2 

which, after a rearrangement of terms and change of variables back to the r — t space, 
becomes the "spherical KdV" equation: 

dpi , dpi 2 dpi ( gy'^po \d^pi pi 

13 



for which a suitable initial condition may be: 

3(u — c 



Pi{r, to) 



-(3-c 



2\-l 



sech 



my 



[U - Cs)CsM 



(r — ut, 



oj 



(56) 



This Gaussian-looking form is motivated by the anal ytic al solution of the KdV equation in 
one dimensional Cartesian coordinates discussed in 
parameters without special meaning. 



15| . Here, the numbers m, tor-- cire 



B. Finite temperature 

Apart from the trivial replacement of (12 7p by flTB]) there is another change when we 
consider nuclear matter at finite temperature. We do not restrict ourselves to the case 
where nuclear matter is saturated. Instead, we shall consider the case where there is an 
equilibrated background with constant density and zero velocity, upon which perturbations 
propagate, but no saturation. The difference is that with saturation, a system is bound and 
more stable, whereas in the present case stability is not guaranteed and this system might 
expand or shrink. In such a situation, perturbations would propagate in an expanding 
medium and the reference density po in dSSD might change with time. This is the scenario 
that we have in heavy ion collisions at RHIC, which we plan to address in the future. Here 
we consider the simpler case of constant po- 

Substituting f|T6|) into f|T3|) we obtain: 

h = ^^PB + d^V'pB (57) 
mv my 

Using the definition of the operators in spherical coordinates we arrive at: 

dh _ gy^ dps ^ ^ gy'^ d^pB ^ ^9v'^ ^^PB ^ 9v'^ 2 dps 

dr my^ dr my^ dr^ my^ r dr"^ my^ dr 

and also 

dh _ gy^ dpB ^ ^9v^ ^ ( ^^Pb \ ^ ^Sv^ ^ ( ^Pb\ ^ 2vdpB 

dt my"^ dt my^ dt \ dr"^ J niy^ r dt\ dr J my^ dr 

Substituting fl58|) and fl59|) into (|TT1) and repeating the steps described in the last section, 
i.e., introducing dimensionless variables, changing variables to the ^ — r space, expanding p 
and V and collecting the terms proportional to a and to o"^ we obtain the following relations 
from the Euler equation: 



^ Ts\ dvi ^ gy^ Po dpi 
Po J niy^ 



(60) 
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and 



a 



/is + 



Ts 



Po 



dv2 dvi dvi 



, gv"^ po dpi 
H ^— Pi— + 



_^ 9v^ Po dp2 _ 
Cs"^ dC, 



dvi gv'^ dpi gv 

at, my^ 



Po d^pi 



After some manipulations, the continuity equation ([7]) is written as: 



dpB ^ dv ^ ^dpB ^ 2pB^ \ ^ ^, 
dt dr dr r J 







(61) 



(62) 



which, after the change of variables and expansion yields the following relations: 

a| --^ + —1=0 

and 



dpi dvi^ 
d^ 



a 



dv2 , dpi 2 dvi dp2 , dpi dvi 2 







(63) 



(64) 



where in the last term of the above expression we made again the approximation ( 150|) . From 
( 160|) and (|63|) we get the relations: 



Ts 

Pb^ 

Po 



9v Po 



and 



Vi = pi 

Substituting these expressions in (I6T1) and (1641) and combining the resulting equations we 
arrive at the finite temperature spherical KdV equation: 

dpi 



dpi ^ ^ dpi 
dt ^ dr 



PB'mv'^Cs^ 



CsPv 



d 



2my2 / dr^ 



d^Pi ^ Pi 

t 







(65) 



25'y^Po y*'"^ dr \z.'iiLv- 
In the numerical studies of this equation we have used an initial condition with the form 
given by with several choices for the parameters. 



One dimensional Cartesian coordinates 



One dimensional perturbations propagating in cold nuclear matter have been discussed 
in detail in 15| and 16| and thus we only give here the obtained differential equation: 

dpi 



upi , ^,/q_ 2^ . dpi ^ ^ " ^ 
dt ^""'dx^ ^""'P' dx + \2Mmv'c 



( ^9v Po ^ d'^^Pi 
dx^ 







(66) 
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which has the following analytical solution: 



my^ {u — Cs)CsM 



9v 



2d Pa 



(x — ut) 



(67) 



At finite temperature we replace fl27|) by (fT6l) and use it directly in (fT3!) without imposing 
any saturation condition. Everything else is the same as described in the last sections. We 
arrive at the following equation: 

dpi ( dcs \ d^pi 



dpi ^ ^ dpi 
dt ^ dx 



2 - c/ - 



2 PBfnv'^Cs'^ \ „ upi 
Cspl- 



dx 



with analytical solution given by: 
3(u — c 



pi{x,t) 



'^9v Po 



2mv^ / dx^ 







{x — ut) 



(68) 



(69) 



V. NUMERICAL RESULTS 



The numerical solution of non-linear differential equations is not very difficult, but may 
be very tricky. We benefited from the reading and hints contained in the textbook |28| . 
which has a special section dedicated to solitons. 

We start our numerical analysis showing in Fig. [T]the solution of the linear KdV equation 
at T = 0, Eq. (1661) . In the upper pannel, Fig. [TK), we use the analytical solution Eq. (1671) 
as initial condition. As expected this pulse propagates without dissipation nor dispersion. 
Any change in the initial condition has noticeable consequences. In Fig. [TJo), we follow the 
evolution of the numerical solution of (166|) for an initial pulse given by (!67|) multiplied by 
a factor two. As it can be seen, the amplitude grows, the width decreases and a second 
bumps appears propagating behind the first. In Fig. [Tfc) we start the evolution with (1671) 
multiplied by a factor seven. Now we have three peaks instead of two. In practical cases, the 
initial conditions will never be exactly those needed to generate a single pulse. Therefore, 
in general we expect to see multiple bumps. In Fig. [1] we can also see that perturbations 
with higher amplitudes propagate faster. 

In Fig. [21 we show the equivalent plot for the spherical case. In contrast to the linear 
case there is a strong damping of the pulse. The dependence on the initial conditions is also 
strong. The main peak very rapidly looses height and develops secondary bumps. 

In Fig. [3] we show for the linear case and for the "optimal" initial condition (167]) the 
evolution of the pulse with time for different temperatures. We can see that, increasing the 
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temperature the pulses move faster and go farther. The same feature can be observed in the 
spherical shown in Fig. |H 

Setting d = in fll4p we recover the standard non-linear Walecka model. In this medium 
the propagation of density perturbations will be governed by the differential equations ([SSD, 
(I65p . (IMj) and (1681) without the terms with the d factor. The absence of the third order 
derivative terms leads to a lack of stability of the solution. The corresponding differential 
equations are "shock wave equations". Out of smooth initial perturbations these equations 
create shock waves. We can see this process in one dimensional Cartesian coordinates in 
Fig. El 

Figure 5a) shows the solution of Eq. (IMj) with d = for the initial condition (1561) . In 
Figs. 5b) and 5c) the initial profile has been multiplied by a factor 2 and 10 respectively. 
In all cases we observe a steepening of the profile until the formation of the shock, followed 
by the dispersion of the wave. We see that the higher is the initial amplitude, the sooner 
the wave breaking and dispersion occur. The solutions of ( !55l) with d = for three different 
initial profiles can be seen in Fig. El where the three pannels, 6a), 6b) and 6c), show the 
evolution of the original profile and then the evolution of this profile multiplied by 2 and by 
10 respectively. In the spherical case the damping is so strong that wave breaking hardly 
happens. 

In Figure [7] we fix one initial profile and study its time evolution for three different 
temperatures. Figs. 7a), 7b) and 7c) show the development of a shock wave at T = 20 
MeV, 70 MeV and 120 MeV respectively. As it can be seen, with increasing temperatures 
the pulse moves faster and the shock formation and the subsequent dispersive breaking occur 
sooner. Fig. El shows the analogous plot for spherical coordinates. 



VI. CONCLUSIONS 



From the equations of relativistic hydrodynamics and with an equation of state obtained 
from a variant of the non-linear Walecka model we have derived a spherical KdV-like equation 
for perturbations in the baryon density. The coefficients of the differential equation are 
determined by the microscopic meson exchange dynamics. This is an improvement over our 
previous study jl8]. Moreover we have included temperature effects and solved numerically 
the resulting differential equations. We have also, for the first time, obtained numerical 
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solutions for the one dimensional Cartesian case at zero and finite temperature. 

The results give a quantitative measure of the dependence of the numerical solutions on 
the initial conditions. We found that, as expected in non-linear problems, the behavior of 
the solutions depends very strongly on the initial conditions. 

The results presented in Fig. [3] are a first step towards a realistic study of the propagation 
of a fast leading particle (coming from a jet) crossing hot hadronic matter. They suggest 
that perturbation may propagate for a relatively long distance preserving features of the 
initial peak structure. This is even more true at higher temperatures. In contrast, in 
the spherical case shown in Fig. HI our results show a strong attenuation, indicating that 
localized perturbations will not survive for long distances. They will instead release energy 
to the medium in a more homogeneous way. This behavior may have consequences for 
astrophysical phenomena and we plan to address this subject in the near future. Switching 
off the cubic derivative term in the Lagrangian density and recovering the standard non- 
linear Walecka model, the propagation of initial density pulses generates shock waves, which 
go through a dispersive breaking. Both the propagation and breaking depend strongly on 
the properties (heigth and width) of the initial pulses and on the temperature of the medium. 
Higher pulses move faster and break earlier. The same effect is observed when we increase 
the temperature. In contrast, spherical pulses are very insensitive to the initial conditions 
and to the temperature. 

We plan to investigate the consequences of our findings both in the relativistic heavy ion 
physics and dense stars physics scenarios. 
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FIG. 1: Time evolution of a density perturbation in arbitrary units in one dimensional Cartesian 
coordinates. The upper pannel shows the evolution of the analytic solution of the KdV equation. 
In the lower pannels we show the evolution of the analytic solution multiplied by a factor 2 and 7 
respectively. 
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FIG. 2: The same as Figure 1 for spherical coordinates. 
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FIG. 3: Time evolution of density perturbations in one dimensional Cartesian coordinates as a 
function of the temperature. The pannels show calculations with temperatures T = 20, 70 and 120 
MeV respectively. 
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FIG. 4: The same as Figure 3 for one dimensional spherical coordinates 
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FIG. 5: Shock wave formation in one dimensional Cartesian coordinates, a) Initial profile with 
heigth 0.25 (arbitrary units); b) with heigth 0.50; c) with heigth 2.5. 
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FIG. 7: Shock wave formation in one dimensional Cartesian coordinates for different temperatures. 
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FIG. 8: The saem as Figure 7 for one dimensional spherical coordinates. 
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